clear 
set more off


********************************************************************************
**********************         USER INPUT        *******************************
********************************************************************************


************ PATHS *************************************************************

local inputdir "${temp}"
local output "${temp}"

**************** Set number of  simulations ***************************

local S = 50000

***********************************************************************
*** OPTIONS BELOW USED TO SELECT WHICH TRANSITION MATRICES TO USE *****

********** SPECIFY depths *********************************************

	local depths 1
	 
********** SPECIFY age bins for moving average *********************************

	local age_bin=2 
	 
************ Set Discount Factor ***********************************************

** Sets discount factor used to generate NPVs  
	** IMPORTANT: Also used within percentile calculations for depth > 1

	local beta  97
	
** options for kernel are "uniform" or "triangular"
	
	local kernel "triangular"

*************** Set starting and ending ages **************************

	** (maxage-minage) must be a multiple of all depths used

local minage = 20

local maxage = 70


*************** Set Type *******************************************************

local types "singleyear"

************ setting list of income variables to use ***************************

local incvars  ptotinc

************ CHOOSE COMPARISON OCCUPATIONS *************************************
	
** Which occupations? ex. - gl ocnums 4 6 - for chiropractors and pharmacists**
	*note : if `s' = "nppes" occnums cannot include 8, 9, or 10, which are ACS only.

local occnums 1 2 3

/*
Dictionary:
1 "Physicians" 

2 "Lawyers" 

*/

************* Set incbin ids **************************************************

	*Should be same as labels used in transition matrices
		* First bin name must be for missing/non-filing
	
local incbins "00" "01" "10" "20" "30" "40" "50" "60" "70" "80" "90" "95" "99"

	* list tops of bins excluding bottom (not filing)
local bintops   .1, .2, .3, .4, .5, .6, .7, .8, .9, .95, .99, 1 

********************************************************************************
#delimit ;

local occs 

"
physicians 
pcps
lawyers 
"
;

#delimit cr


foreach depvar in `incvars'{
	foreach i in `occnums'{

	local o: word `i' of `occs'
		foreach d in `depths'{
			foreach f in "`types'"{

				
				cap use "`inputdir'/tm_`o'_`depvar'_`d'years_`f'_`kernel'.dta", clear

				if _rc == 0 {
					
					*recoding income bins
					gen current_bin = .
					gen future_bin = . 
					local bin_num = -1
					foreach bin in "`incbins'" {
						local bin_num = `bin_num' + 1
						replace current_bin = `bin_num' if initial == "`bin'"
						replace future_bin = `bin_num' if subsequent == "`bin'"
					}

					sort age current_bin future_bin

					gen obsnumber = _n
					tsset obsnumber

					gen cum_prob = .
					replace cum_prob = tr_prob in 1

					local maxobs = _N
					
					*generating conditional CMFs 
					forvalues n = 2(1)`maxobs'{
						qui replace cum_prob = L.cum_prob + tr_prob in `n'
						qui replace cum_prob = tr_prob if future_bin == 0
					}

					*number of simulations, set  in user inputs
					set obs `S'
					*Simulating Trajectories
					forvalues a = `minage'(`d')`maxage'{
						local a_1 = `a' - `d'
						* Have to generate starting distribution at youngest age
						if `a' == `minage'{
							gen bin_`a' = .
							*Determining proportion of non-filers at youngest age category
							qui total tr_count if age == `a'  & current_bin== 0
							matrix tot_notfiling = e(b)
							local tot_notfiling = tot_notfiling[1,1]
							qui total tr_count if age == `a'
							matrix tot = e(b)
							local tot = tot[1,1]
							
							local howmany = floor(`S'*`tot_notfiling'/`tot')
							if `howmany'==0 local howmany = 1
							
							* filling in nonfilers
							replace bin_`a' = -1 in 1/`howmany'
							
							
							* categorizing filers
							gen draw = (_n - `howmany')/(`S' - `howmany') if bin_`a' != -1
							replace bin_`a' = irecode(draw , `bintops' ) if bin_`a' != -1
							replace bin_`a' = bin_`a' + 1
							drop draw
						}
						* Simulating the strings of subsequent income bins.
						
						else{
							gen draw = runiform()
							gen bin_`a' = .

							forvalues prev_bin = 0(1)`bin_num'{

									
									qui levelsof cum_prob if current_bin == `prev_bin' & age == `a_1', local(cuts) sep(,)
									*display `cuts'
									replace bin_`a' = irecode(draw, `cuts') if bin_`a_1' == `prev_bin'
										
								

							}
							drop draw

						}

						
					}
					
					*Calculating NPVs
					foreach b in `beta'{
					
					gen NPV_d0`b' = 0
					*equivalent per-period annuity income (abstracting from uncertainty)
					gen const_equiv_d0`b' = .
					
					forvalues a = `minage'(`d')`maxage'{
						gen inc_`a'_`b' = 0
						*extracting mean income by initial bin. 
						forvalues bin = 0(1)`bin_num'{
							*adding in incomes for intermediate ages (a_int) if `d' > 1 
							local a_bintop = `a'+`d' - 1
							forvalues a_int = `a'/`a_bintop'{
								*Only do this for intermediate years if "singleyear" type, since multiyears are filled in within transition matrix code. 
								*Also limiting so don't run over maximum age (Otherwise will require alteration to definition of `T', below).
								if ((`a_int' == `a' | "`f'" == "singleyear") & `a_int' <= `maxage'){
									qui mean mean_`depvar' if age == `a_int' & current_bin == `bin' [aw = tr_prob]
									matrix eb = e(b)
									local meaninc = eb[1,1]
									replace inc_`a'_`b' = inc_`a'_`b' + 0.`b'^(`a_int' - `a')*`meaninc' if bin_`a' == `bin'
								}
							}
							
						}
						* updating NPV
						replace NPV_d0`b' = NPV_d0`b' + 0.`b'^(`a'- `minage' ) * inc_`a'_`b'
						
						disp "`b'" "`a'" "`d'" "`f'" "`o'" "`depvar'"
						
					}
					local T = `maxage' - `minage'
					replace const_equiv_d0`b' = NPV_d0`b'/((1 - 0.`b'^`T')/(1 -0.`b'))
					
					
					preserve
					keep occupation incvar depth f bin_* inc_* NPV_d0`b' const_equiv_d0`b'
					
					
					*replace incvar = "`depvar'"
					replace occupation = "`o'"
					replace depth = `d'
					replace f = "`f'"
					
					save "`output'/simulated_npvs_`o'_`depvar'_`d'years_`f'_`kernel'_b`b'.dta", replace
					restore
					
					}
				}
			}
		}
	}
}

